#------------------------------------- #
# The Journal of Politics
# When Do Citizens Resist The Use of AI Algorithms in Public Policy?
#
# Replication Code for Decision Maker Experiment (appears in the manuscript under the title: Policy Evaluation Experiment))
# 
# Author: Shir Raviv
# Date: Feb 19, 2025
#------------------------------------- #
# Clear environment and free memory ---------------
rm(list = ls())
gc()

# Load required packages -----------
required_packages <- c("dplyr", "readr", 'arsenal','tidyr','estimatr','patchwork',
                       "expss","scales",'ggeffects','lme4','forcats',
                       "vtable","xtable",'stargazer',"lme4",'texreg','broom','rlang',
                       "rio", "broom","ggplot2")
for (pkg in required_packages) {
  if (!require(pkg, character.only = TRUE)) {
    install.packages(pkg, dependencies = TRUE)
    library(pkg, character.only = TRUE)
  }
}


# Set working directory and define export paths ---------------
#setwd("insert here.. /JOP Replication")
plotpath <- "figures/"

# Load Survey Data  ---------------
df <- readRDS("DM_wide_data.RDS")
# Table (A-15) Balance Tables----

# Define the list of demographic variables
demographics <- c("female_fac", "age_cat", "low_edu_somecol_cat", "white","pid3_fac", "high_diglit")



# Balance table for Homelessness

balance_homless <- df %>%
  filter(DMexp_order_homless_first == 1) %>%
  filter(DMexp_treatment_homless != "ADS+HDS") %>%
  select(all_of(demographics), DMexp_treatment_homless)

balance_homless_tab <- CreateTableOne(vars = demographics,
                                      strata = "DMexp_treatment_homless",
                                      data = balance_homless,
                                      factorVars = demographics)
print(balance_homless_tab, showAllLevels = TRUE)



# Balance table for education

balance_educ <- df %>%
  filter(DMexp_order_educ_first == 1) %>%
  filter(DMexp_treatment_educ != "ADS+HDS") %>%
  select(all_of(demographics), DMexp_treatment_educ)

balance_educ_tab <- CreateTableOne(vars = demographics,
                                   strata = "DMexp_treatment_educ",
                                   data = balance_educ,
                                   factorVars = demographics)
print(balance_educ_tab, showAllLevels = TRUE)



# Balance table for policing

balance_police <- df %>%
  filter(DMexp_order_polic_first == 1) %>%
  filter(DMexp_treatment_polic != "ADS+HDS") %>%
  select(all_of(demographics), DMexp_treatment_polic)

balance_police_tab <- CreateTableOne(vars = demographics,
                                     strata = "DMexp_treatment_polic",
                                     data = balance_police,
                                     factorVars = demographics)
print(balance_police_tab, showAllLevels = TRUE)



# Balance table for child welfare

balance_child <- df %>%
  filter(DMexp_order_child_first == 1) %>%
  filter(DMexp_treatment_child != "ADS+HDS") %>%
  select(all_of(demographics), DMexp_treatment_child)

balance_child_tab <- CreateTableOne(vars = demographics,
                                    strata = "DMexp_treatment_child",
                                    data = balance_child,
                                    factorVars = demographics)
print(balance_child_tab, showAllLevels = TRUE)



# Figure (5) Average policy support, by decision-maker and context ----------

#  Data prep for top panel (policy support) 
mean_DMexp_bin_outcome_homless<-
  df %>%   
  filter(DMexp_order_homless_first==1)%>%   
  dplyr::select(DMexp_treatment_homless, DMexp_bin_outcome_homless) %>%
  filter(!is.na(DMexp_bin_outcome_homless)) %>%
  filter(!is.na(DMexp_treatment_homless)) %>%
  group_by(DMexp_treatment_homless) %>%     
  summarize(mean=mean(DMexp_bin_outcome_homless), n=n(), sd=sd(DMexp_bin_outcome_homless), se=sd/sqrt(n)) %>%   
  mutate(DM=DMexp_treatment_homless)

mean_DMexp_bin_outcome_child<-
  df %>%  
  filter(DMexp_order_child_first==1)%>%   
  dplyr::select(DMexp_treatment_child, DMexp_bin_outcome_child) %>%
  filter(!is.na(DMexp_treatment_child)) %>%
  filter(!is.na(DMexp_bin_outcome_child)) %>%
  group_by(DMexp_treatment_child) %>%     
  summarize(mean=mean(DMexp_bin_outcome_child), n=n(), sd=sd(DMexp_bin_outcome_child), se=sd/sqrt(n))%>%   
  mutate(DM=DMexp_treatment_child)     
mean_DMexp_bin_outcome_educ<-
  df %>%   
  filter(DMexp_order_educ_first==1 )%>%   
  dplyr::select(DMexp_treatment_educ, DMexp_bin_outcome_educ) %>%
  filter(!is.na(DMexp_treatment_educ)) %>%
  filter(!is.na(DMexp_bin_outcome_educ)) %>%
  group_by(DMexp_treatment_educ) %>%     
  summarize(mean=mean(DMexp_bin_outcome_educ), n=n(), sd=sd(DMexp_bin_outcome_educ), se=sd/sqrt(n)) %>%   
  mutate(DM=DMexp_treatment_educ)    


mean_DMexp_bin_outcome_polic<-
  df %>%   
  filter(DMexp_order_polic_first==1)%>%   
  dplyr::select(DMexp_treatment_polic, DMexp_bin_outcome_polic) %>%
  filter(!is.na(DMexp_treatment_polic)) %>%
  filter(!is.na(DMexp_bin_outcome_polic)) %>%
  group_by(DMexp_treatment_polic) %>%     
  summarize(mean=mean(DMexp_bin_outcome_polic), n=n(), sd=sd(DMexp_bin_outcome_polic), se=sd/sqrt(n)) %>%   
  mutate(DM=DMexp_treatment_polic)  

df_fig5 <- bind_rows(
  mean_DMexp_bin_outcome_homless,
  mean_DMexp_bin_outcome_educ,
  mean_DMexp_bin_outcome_child,
  mean_DMexp_bin_outcome_polic,
  .id="context"
)

df_fig5$context <- factor(df_fig5$context,
                          levels = c("1", "2", "3", "4"),  
                          labels = c("Public Housing", 
                                     "Funding Educational Programs", 
                                     "Child Abuse Allegations", 
                                     "Police Patrols")
)

df_fig5$DM <- factor(df_fig5$DM,
                     levels = c("ADS","HDS"),  
                     labels = c("Algorithmic\nDecision-Maker","Human\nDecision-Maker")
)


fig5_bars <- df_fig5 %>%
  filter(DM != "ADS+HDS") %>%  
  ggplot(aes(x = DM, y = mean, fill = DM)) +
  geom_col(position = position_dodge(width=0.6), width = 0.5) +
  geom_errorbar(aes(ymin = mean - 1.96*se, ymax = mean + 1.96*se),
                width = 0, color="black", size = 0.5,
                position = position_dodge(width=0.6)) +
  scale_y_continuous(labels = percent_format(accuracy = 1),
                     limits = c(0, 0.70)) +  # or adjust as needed
  scale_fill_manual(values = c("#406060","#99A8A3")) + 
  facet_wrap(~ context, ncol=2) +
  labs(x = NULL, y = "Support") +
  theme_minimal(base_family = "Times New Roman", base_size = 14) +
  theme(
    panel.background = element_rect(fill = "white", colour = "gray40"),
    legend.position = "none",
    strip.text = element_text(size = 12),
    axis.title.y = element_text(vjust = 1.5, size=11.5)
  )
fig5_bars
# Data prep for bottom panel (difference in means)
homless_HDM_ADM <- difference_in_means(formula = DMexp_bin_outcome_homless ~ 
                                         DMexp_treatment_homless, 
                                       data = subset(df,DMexp_order_homless_first==1 ), 
                                       condition1 = "HDS", 
                                       condition2 = "ADS")

child_HDM_ADM <- difference_in_means(formula = DMexp_bin_outcome_child ~ 
                                       DMexp_treatment_child, 
                                     data = subset(df,DMexp_order_child_first==1 ), 
                                     condition1 = "HDS", 
                                     condition2 = "ADS")

educ_HDM_ADM <- difference_in_means(formula = DMexp_bin_outcome_educ ~ 
                                      DMexp_treatment_educ, 
                                    data = subset(df,DMexp_order_educ_first==1 ), 
                                    condition1 = "HDS", 
                                    condition2 = "ADS")

police_HDM_ADM <- difference_in_means(formula = DMexp_bin_outcome_polic ~ 
                                        DMexp_treatment_polic, 
                                      data = subset(df, DMexp_order_polic_first==1 ), 
                                      condition1 = "HDS", 
                                      condition2 = "ADS")

m.police_HDM_ADM <-  tidy(police_HDM_ADM, conf.int = TRUE)
m.educ_HDM_ADM <-  tidy(educ_HDM_ADM, conf.int = TRUE)
m.child_HDM_ADM <-  tidy(child_HDM_ADM, conf.int = TRUE)
m.homless_HDM_ADM <-  tidy(homless_HDM_ADM, conf.int = TRUE)



fig5_df_coef <- bind_rows(
  m.homless_HDM_ADM,
  m.educ_HDM_ADM,
  m.police_HDM_ADM,
  m.child_HDM_ADM,
  .id="domain"
) %>%
  select(-c(term, outcome))

fig5_df_coef$domain <- factor(fig5_df_coef$domain, 
                              levels=c("1","2","3","4"),
                              labels=c("Public Housing", "Education", "Policing", "Child Welfare")
)

# Plot the bottom panel (coefficient plot)
fig5_df_coef$domain <- factor(
  fig5_df_coef$domain,
  levels = c("Public Housing", "Child Welfare", "Education", "Policing")
)

fig5_coef <- fig5_df_coef %>%
  ggplot(aes(y = domain, x = estimate)) +
  geom_vline(xintercept = 0, linetype = 2, color = "#ACBBBF", size = 1.2) +
  geom_point(size = 2.5, color = "#406060") +
  geom_errorbarh(aes(xmin = estimate - 1.96 * std.error, 
                     xmax = estimate + 1.96 * std.error),
                 height = 0, 
                 color = "#406060", 
                 alpha = 0.7, 
                 size = 1) +
  geom_text(aes(label = sprintf("%.3f", estimate)),
            nudge_y = 0.25,     # shift the text slightly
            size = 3, 
            color = "black") + coord_flip()+
  scale_x_continuous(limits = c(-0.27, 0.27)) +
  theme_minimal(base_family = "Times New Roman", base_size = 14) +
  labs(x = "", y="") +
  theme(
    axis.text.y = element_blank(),
    axis.ticks.x= element_blank(),
    legend.position = "none",
    panel.background = element_rect(fill = "white", colour = "gray40"),
    plot.background = element_rect(fill = "white", colour = "white"),
    axis.title.y = element_blank(),  # no label on y-axis
    axis.title.x = element_text(size = 12) ) 

fig5_coef

# Combine the two panels (top & bottom) into one figure:
fig5 <- fig5_bars / fig5_coef + 
  plot_layout(heights = c(3, 1))  # top panel taller than bottom

fig5


ggsave(fig5,
       filename = paste0(plotpath, "fig5.png"),
       device = "png",
       height = 7.5, width = 7, dpi = 600)


ggsave(fig5,
       filename = paste0(plotpath, "fig5.tif"),
       device = "tiff",
       height = 7.5, width = 7, dpi = 600)


# Table (A-16) Effects of ADS on the Support for Policy Proposals -------
kable(as.data.frame(fig5_df_coef), digits=3, "latex", booktabs=TRUE)

# Table (A-17) ADS versus HDM assisted by ADS ------
homless_Hybrid <- difference_in_means(formula = DMexp_bin_outcome_homless ~ 
                                            DMexp_treatment_homless, 
                                          data = subset(df,DMexp_order_homless_first==1 ), 
                                          condition2 = "ADS",
                                          condition1 = "ADS+HDS")

child_Hybrid <- difference_in_means(formula = DMexp_bin_outcome_child ~ 
                                          DMexp_treatment_child, 
                                        data = subset(df,DMexp_order_child_first==1 ), 
                                        condition2 = "ADS",
                                        condition1 = "ADS+HDS")

educ_Hybrid <- difference_in_means(formula = DMexp_bin_outcome_educ ~ 
                                         DMexp_treatment_educ, 
                                       data = subset(df,DMexp_order_educ_first==1 ), 
                                       condition2 = "ADS",
                                       condition1 = "ADS+HDS")

police_Hybrid <- difference_in_means(formula = DMexp_bin_outcome_polic ~ 
                                           DMexp_treatment_polic, 
                                         data = subset(df, DMexp_order_polic_first==1 ), 
                                         condition2 = "ADS",
                                         condition1 = "ADS+HDS")

m.homless_Hybrid <-  tidy(homless_Hybrid, conf.int = TRUE)
m.educ_Hybrid <-  tidy(educ_Hybrid, conf.int = TRUE)
m.police_Hybrid <-  tidy(police_Hybrid, conf.int = TRUE)
m.child_Hybrid <-  tidy(child_Hybrid, conf.int = TRUE)


df_coef_data_Hybrid <- bind_rows(
  m.homless_Hybrid, 
  m.educ_Hybrid,
  m.police_Hybrid,
  m.child_Hybrid,
 .id="domain")

df_coef_data_Hybrid$domain<-factor(df_coef_data_Hybrid$domain, 
                                         levels=c(1,2,3,4), 
                                         labels=c( "Public Housing", "Education",
                                                   "Policing", "Child Welfare"))

df_coef_data_Hybrid<-df_coef_data_Hybrid%>%
  dplyr::select(-c(term, outcome))

df_coef_data_Hybrid$domain<-factor(df_coef_data_Hybrid$domain)
kable(as.data.frame(df_coef_data_Hybrid), digits=3, "latex", booktabs=TRUE)



# Table (A-18) Additional Results -------

# Create factor variables for treatments and set "HDS" as the reference level.

df$DMexp_treatment_homless2<-as.factor(df$DMexp_treatment_homless)
df$DMexp_treatment_homless2<-relevel(df$DMexp_treatment_homless2, 
                                     ref="HDS")
df$DMexp_treatment_child2<-as.factor(df$DMexp_treatment_child)
df$DMexp_treatment_child2<-relevel(df$DMexp_treatment_child2, 
                                   ref="HDS")
df$DMexp_treatment_educ2<-as.factor(df$DMexp_treatment_educ)
df$DMexp_treatment_educ2<-relevel(df$DMexp_treatment_educ2, 
                                  ref="HDS")
df$DMexp_treatment_polic2<-as.factor(df$DMexp_treatment_polic)
df$DMexp_treatment_polic2<-relevel(df$DMexp_treatment_polic2, 
                                   ref="HDS")

# -- Columns 1-4 in Table A-18 (Public Housing domain)  -----
#   - Model 1: Basic treatment effect
lm1_housing<-lm(DMexp_bin_outcome_homless~ DMexp_treatment_homless2,
                data=subset(df, DMexp_order_homless_first==1))
#   - Model 2: Treatment effect with demographic and attitudinal covariates
lm2_housing<-lm(DMexp_bin_outcome_homless~ DMexp_treatment_homless2+
                  female_fac+birthyear+low_edu_somecol_cat+white+pid3_fac+
                  high_diglit+inattentive,
                data=subset(df, DMexp_order_homless_first==1))
#   - Model 3: Treatment effect with treatment order
lm3_housing<-lm(DMexp_bin_outcome_homless~ DMexp_treatment_homless2+DMexp_order_homless,
                data=subset(df))
#   - Model 4: Treatment effect with order and covariates
lm4_housing<-lm(DMexp_bin_outcome_homless~ DMexp_treatment_homless2+DMexp_order_homless+
                  female_fac+birthyear+low_edu_somecol_cat+white+pid3_fac+
                  high_diglit+inattentive,
                data=subset(df))
# Columns 1-4 of Table A-18 (Public Housing domain) 
stargazer(lm1_housing, lm2_housing, lm3_housing, lm4_housing, 
          type="text", 
          dep.var.labels = "Support (public housing)",
          covariate.labels=c('ADS', 'HDM assisted by ADS', 
                             'Female', 'Age','Some college or less', 'White',
                             'Independent', 'Republican', 
                             'High Tech Literacy', "Inattentive", 'Order', 'Constant'),
          omit.stat = c('adj.rsq', 'ser', 'f'),
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001), star.char = c('†', '*', '**', '***'),
          notes = '†p<0.1; *p<0.05; **p<0.01; ***p<0.001', notes.append = F)

# -- Columns 5-8 in Table A-18 (Child Welfare domain) --------
#   - Model 1: Basic treatment effect
lm1_childwel<-lm(DMexp_bin_outcome_child~ DMexp_treatment_child2,
                 data=subset(df, DMexp_order_child_first==1))
#   - Model 2: Treatment effect with demographic and attitudinal covariates
lm2_childwel<-lm(DMexp_bin_outcome_child~ DMexp_treatment_child2+
                   female_fac+birthyear+low_edu_somecol_cat+white+pid3_fac+
                   high_diglit+inattentive,
                 data=subset(df, DMexp_order_child_first==1))
#   - Model 3: Treatment effect with treatment order
lm3_childwel<-lm(DMexp_bin_outcome_child~ DMexp_treatment_child2+DMexp_order_child,
                 data=subset(df))
#   - Model 4: Treatment effect with order and covariates
lm4_childwel<-lm(DMexp_bin_outcome_child~ DMexp_treatment_child2+DMexp_order_child+
                   female_fac+birthyear+low_edu_somecol_cat+white+pid3_fac+
                   high_diglit+inattentive,
                 data=subset(df))
# Columns 5-8 of Table A-18 
stargazer(lm1_childwel, lm2_childwel, lm3_childwel, lm4_childwel, 
          type="text", 
          dep.var.labels = "Support (Child Welfare)",
          covariate.labels=c('ADS', 'HDM assisted by ADS', 
                             'Female', 'Age','Some college or less', 'White',
                             'Independent', 'Republican', 
                             'High Tech Literacy', "Inattentive", 'Order', 'Constant'),
          omit.stat = c('adj.rsq', 'ser', 'f'),
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001), star.char = c('†', '*', '**', '***'),
          notes = '†p<0.1; *p<0.05; **p<0.01; ***p<0.001', notes.append = F)



# -- Columns 9-12 in Table A-18 (Education domain) --------
#   - Model 1: Basic treatment effect
lm1_educ<-lm(DMexp_bin_outcome_educ~ DMexp_treatment_educ2,
             data=subset(df, DMexp_order_educ_first==1))
#   - Model 2: Treatment effect with demographic and attitudinal covariates
lm2_educ<-lm(DMexp_bin_outcome_educ~ DMexp_treatment_educ2+
               female_fac+birthyear+low_edu_somecol_cat+white+pid3_fac+
               high_diglit+inattentive,
             data=subset(df, DMexp_order_educ_first==1))
#   - Model 3: Treatment effect with treatment order
lm3_educ<-lm(DMexp_bin_outcome_educ~ DMexp_treatment_educ2+DMexp_order_edu,
             data=subset(df))
#   - Model 4: Treatment effect with order and covariates
lm4_educ<-lm(DMexp_bin_outcome_educ~ DMexp_treatment_educ2+DMexp_order_edu+
               female_fac+birthyear+low_edu_somecol_cat+white+pid3_fac+
               high_diglit+inattentive,
             data=subset(df))
# Columns 9-12 of Table A-18 
stargazer(lm1_educ, lm2_educ, lm3_educ, lm4_educ, 
          type="text", 
          dep.var.labels = "Support (Education)",
          covariate.labels=c('ADS', 'HDM assisted by ADS', 
                             'Female', 'Age','Some college or less', 'White',
                             'Independent', 'Republican', 
                             'High Tech Literacy', "Inattentive", 'Order', 'Constant'),
          omit.stat = c('adj.rsq', 'ser', 'f'),
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001), star.char = c('†', '*', '**', '***'),
          notes = '†p<0.1; *p<0.05; **p<0.01; ***p<0.001', notes.append = F)

# -- Columns 13-16 in Table A-18 (Policing domain) -----
#   - Model 1: Basic treatment effect
lm1_polic<-lm(DMexp_bin_outcome_polic~ DMexp_treatment_polic2,
              data=subset(df, DMexp_order_polic_first==1))
#   - Model 2: Treatment effect with demographic and attitudinal covariates
lm2_polic<-lm(DMexp_bin_outcome_polic~ DMexp_treatment_polic2+
                female_fac+birthyear+low_edu_somecol_cat+white+pid3_fac+
                high_diglit+inattentive,
              data=subset(df, DMexp_order_polic_first==1))
#   - Model 3: Treatment effect with treatment order
lm3_polic<-lm(DMexp_bin_outcome_polic~ DMexp_treatment_polic2+DMexp_order_polic,
              data=subset(df))
#   - Model 4: Treatment effect with order and covariates
lm4_polic<-lm(DMexp_bin_outcome_polic~ DMexp_treatment_polic2+DMexp_order_polic+
                female_fac+birthyear+low_edu_somecol_cat+white+pid3_fac+
                high_diglit+inattentive,
              data=subset(df))
# Columns 13-16 of Table A-18
stargazer(lm1_polic, lm2_polic, lm3_polic, lm4_polic, 
          type="text", 
          dep.var.labels = "Support (Policing)",
          covariate.labels=c('ADS', 'HDM assisted by ADS', 
                             'Female', 'Age','Some college or less', 'White',
                             'Independent', 'Republican', 
                             'High Tech Literacy', "Inattentive", 'Order', 'Constant'),
          omit.stat = c('adj.rsq', 'ser', 'f'),
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001), star.char = c('†', '*', '**', '***'),
          notes = '†p<0.1; *p<0.05; **p<0.01; ***p<0.001', notes.append = F)



# Table (A-19) Alternative Measures of Outcomes --------
# Reverse the 5-point outcome scales so that higher values indicate higher support
df$DMexp_outcome_homless_rev<-6-df$DMexp_outcome_homless
df$DMexp_outcome_child_rev<-6-df$DMexp_outcome_child
df$DMexp_outcome_educ_rev<-6-df$DMexp_outcome_educ
df$DMexp_outcome_polic_rev<-6-df$DMexp_outcome_polic

# Generate binary outcome variables: "Strongly Support" when outcome == 1, else 0.
df <- df %>%
  mutate(DMexp_bin_strong_outcome_homless = case_when(
    DMexp_outcome_homless == 1 ~ 1,
    DMexp_outcome_homless %in% c(2, 3, 4, 5) ~ 0)) %>%
  mutate(DMexp_bin_strong_outcome_child = case_when(
    DMexp_outcome_child == 1 ~ 1,
    DMexp_outcome_child %in% c(2, 3, 4, 5) ~ 0)) %>%
  mutate(DMexp_bin_strong_outcome_educ = case_when(
    DMexp_outcome_educ == 1 ~ 1,
    DMexp_outcome_educ %in% c(2, 3, 4, 5) ~ 0)) %>%
  mutate(DMexp_bin_strong_outcome_polic = case_when(
    DMexp_outcome_polic == 1 ~ 1,
    DMexp_outcome_polic %in% c(2, 3, 4, 5) ~ 0))

#   - A model using the reversed 5-point outcome (continuous)
lm1_housing<-lm(DMexp_outcome_homless_rev~ DMexp_treatment_homless2,
                data=subset(df, DMexp_order_homless_first==1))
lm2_housing<-lm(DMexp_bin_strong_outcome_homless~ DMexp_treatment_homless2,
                data=subset(df, DMexp_order_homless_first==1))

#   - A model using the reversed 5-point outcome (continuous)
lm1_childwel<-lm(DMexp_outcome_child_rev~ DMexp_treatment_child2,
                 data=subset(df, DMexp_order_child_first==1))
#   - A model using the binary "strong support" outcome
lm2_childwel<-lm(DMexp_bin_strong_outcome_child~ DMexp_treatment_child2,
                 data=subset(df, DMexp_order_child_first==1))

#   - A model using the reversed 5-point outcome (continuous)
lm1_educ<-lm(DMexp_outcome_educ_rev~ DMexp_treatment_educ2,
             data=subset(df, DMexp_order_educ_first==1))
#   - A model using the binary "strong support" outcome
lm2_educ<-lm(DMexp_bin_strong_outcome_educ~ DMexp_treatment_educ2,
             data=subset(df, DMexp_order_educ_first==1))

#   - A model using the reversed 5-point outcome (continuous)
lm1_polic<-lm(DMexp_outcome_polic_rev~ DMexp_treatment_polic2,
              data=subset(df, DMexp_order_polic_first==1))
#   - A model using the binary "strong support" outcome
lm2_polic<-lm(DMexp_bin_strong_outcome_polic~ DMexp_treatment_polic2,
              data=subset(df, DMexp_order_polic_first==1))

coef_names_order1 <- c('Constant','ADS', 'HDM assisted by ADS')
names(lm1_housing$coefficients) <- coef_names_order1
names(lm2_housing$coefficients) <- coef_names_order1
names(lm2_childwel$coefficients) <- coef_names_order1
names(lm1_childwel$coefficients) <- coef_names_order1
names(lm2_educ$coefficients) <- coef_names_order1
names(lm1_educ$coefficients) <- coef_names_order1
names(lm2_polic$coefficients) <- coef_names_order1
names(lm1_polic$coefficients) <- coef_names_order1

# LaTeX table
models <- list(
  lm1_housing, lm1_childwel, lm1_educ, lm1_polic,
  lm2_housing, lm2_childwel, lm2_educ, lm2_polic)
print(names(coef(models[[1]])))

texreg(models,
       custom.model.names = c(
         "Housing (5-pt)", "Child (5-pt)", "Educ (5-pt)", "Polic (5-pt)",
         "Housing (bin)", "Child (bin)", "Educ (bin)", "Polic (bin)"),
       stars = c(0.1, 0.05, 0.01, 0.001),digits = 3, 
      
       custom.note = "†p<0.1; *p<0.05; **p<0.01; ***p<0.001")  

# Table (A-20) Interaction - Policy context and Decision maker --------
df<-df%>%
  mutate(order_first=case_when(
    DMexp_order_homless_first==1~"Public Housing",
    DMexp_order_child_first==1~"Child Welfare",
    DMexp_order_educ_first==1~"Education",
    DMexp_order_polic_first==1~"Policing"))
df$order_first<-factor(df$order_first)

df<-df%>%
  mutate(type_collective=case_when(
    order_first=="Public Housing" |
      order_first=="Child Welfare"~ 0,
    order_first=="Education" |
      order_first=="Policing"~ 1))
df$type_collective<-factor(df$type_collective)


df<-df%>%
  mutate(type_sanctioning=case_when(
    order_first=="Public Housing" |
      order_first=="Education" ~ 0,
    order_first=="Child Welfare" |
      order_first=="Policing"~ 1))
df$type_sanctioning<-factor(df$type_sanctioning)


df<-df%>%
  mutate(DM_bin_outcome_first=case_when(
    DMexp_order_homless_first==1~DMexp_bin_outcome_homless,
    DMexp_order_child_first==1~DMexp_bin_outcome_child,
    DMexp_order_educ_first==1~DMexp_bin_outcome_educ,
    DMexp_order_polic_first==1~DMexp_bin_outcome_polic))

df<-df%>%
  mutate(DM_outcome_rev_first=case_when(
    DMexp_order_homless_first==1~DMexp_outcome_homless_rev,
    DMexp_order_child_first==1~DMexp_outcome_child_rev,
    DMexp_order_educ_first==1~DMexp_outcome_educ_rev,
    DMexp_order_polic_first==1~DMexp_outcome_polic_rev))

df<-df%>%
  mutate(dm_treatment=case_when(
    DMexp_order_homless_first==1~DMexp_treatment_homless2,
    DMexp_order_child_first==1~DMexp_treatment_child2,
    DMexp_order_educ_first==1~DMexp_treatment_educ2,
    DMexp_order_polic_first==1~DMexp_treatment_polic2))


lm1<- lm(DM_bin_outcome_first~ dm_treatment+order_first+dm_treatment*order_first, data=df)
lm2<- lm(DM_bin_outcome_first~ dm_treatment+order_first+dm_treatment*order_first+
           birthyear+female_fac+pid3_fac+low_edu_somecol_cat+white+high_diglit, data=df)

lm3<- lm(DM_outcome_rev_first~ dm_treatment+order_first+dm_treatment*order_first, data=df)
lm4<- lm(DM_outcome_rev_first~ dm_treatment+order_first+dm_treatment*order_first+
           birthyear+female_fac+pid3_fac+low_edu_somecol_cat+white+high_diglit, data=df)



stargazer(lm1,lm2, lm3,lm4,
          title="Interaction - Policy context and Decision maker",
          type="latex", 
          dep.var.labels = c("(binary, 2 last cat)","(5-point scale)"),
          order = c(
            "dm_treatmentADS:order_firstEducation",
            "dm_treatmentADS:order_firstPolicing",
            "dm_treatmentADS:order_firstPublic Housing",
            "dm_treatmentADS",
            "order_firstEducation",
            "order_firstPolicing",
            "order_firstPublic Housing",
            "dm_treatmentADS\\+HDS:order_firstEducation",
            "dm_treatmentADS\\+HDS:order_firstPolicing",
            "dm_treatmentADS\\+HDS:order_firstPublic Housing",
            "dm_treatmentADS\\+HDS",
            
            # Demographic variables
            "birthyear",
            "female_facFemale",
            "pid3_facIndependent",
            "pid3_facRepublican",
            "low_edu_somecol_catSome College or Less",
            "whiteWhite",
            "high_diglitHigh Digital Literacy" ),
          covariate.labels = c( 'ADS X Education',   'ADS X Policing',   'ADS X Public Housing','ADS',
                                'Education', 'Policing', 'Public Housing',
                                'ADS+HDS X Education','ADS+HDS X Policing', 'ADS+HDS X Public Housing',
                                'ADS+HDS',
                                'Age','Female','Independent','Republican',
                                'Some college or less','White', 'High Digital Literacy'),
          
          omit.stat = c('adj.rsq', 'ser', 'f'),
          star.cutoffs = c(0.1, 0.05, 0.01, 0.001), star.char = c('†', '*', '**', '***'),
          notes = '†p<0.1; *p<0.05; **p<0.01; ***p<0.001', notes.append = F)





                                           
        